In [1]:
import scvelo as scv
scv.logging.print_version()
Running scvelo 0.2.1 (python 3.7.6) on 2020-06-04 16:18.
In [2]:
scv.settings.verbosity = 3  # show errors(0), warnings(1), info(2), hints(3)
scv.settings.presenter_view = True  # set max width size for presenter view
scv.set_figure_params('scvelo')  # for beautified visualization
  • The analysis is based on the in-built pancreas data.
  • To run velocity analysis on your own data, read your file (loom, h5ad, csv …) to an AnnData object with adata = scv.read('path/file.loom', cache=True).
  • If you want to merge your loom file into an already existing AnnData object, use scv.utils.merge(adata, adata_loom).
In [13]:
adata = scv.datasets.pancreas()
adata
Out[13]:
AnnData object with n_obs × n_vars = 3696 × 27998
    obs: 'clusters_coarse', 'clusters', 'S_score', 'G2M_score'
    var: 'highly_variable_genes'
    uns: 'clusters_coarse_colors', 'clusters_colors', 'day_colors', 'neighbors', 'pca'
    obsm: 'X_pca', 'X_umap'
    layers: 'spliced', 'unspliced'
    obsp: 'distances', 'connectivities'

Names of observations and variables can be accessed via adata.obs_names and adata.var_names, respectively. AnnData objects can be sliced like dataframes, for example, adata_subset = adata[:, list_of_gene_names].

In [11]:
scv.pl.proportions(adata)
  • ここでは、スプライスされた/されていない数の割合が表示されます。使用したプロトコル(Drop-Seq、Smart-Seq)にもよりますが、通常、intron配列を含むスプライスされていない分子が10%~25%程度存在します。
  • また、スプライシング効率の一貫性を確認するために、クラスターレベルでのばらつきを調べることをお勧めします。ここでは、予想通りのばらつきが見られ、周期性管細胞 (cycling ductal cells)ではスプライスされていない分子の割合がわずかに低く、その後、多くの遺伝子が転写され始めるNgn3-high細胞やPre-endocrine細胞では、細胞運命コミットメントでの割合が高くなっていることがわかりました。

Preprocess the Data

  • 遺伝子選択(最小発現量とhigh variability (dispersion))、
  • 細胞をTotal countで正規化し、対数化

フィルタリングと正規化は、スプライスされた/されていないカウントにも同じように適用されます。

これらはすべて単一の関数 scv.pp.filter_and_normalize にまとめられており、以下のコマンドが含まれます。

In [12]:
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
#scv.pp.filter_genes(adata, min_shared_counts=20)
#scv.pp.normalize_per_cell(adata)
#scv.pp.filter_genes_dispersion(adata, n_top_genes=2000)
#scv.pp.log1p(adata)
Filtered out 20801 genes that are detected 20 counts (shared).
Normalized count data: X, spliced, unspliced.
  • さらに、PCA空間の最近傍点間で計算された1次と2次のモーメント(平均と非中心分散)が必要で、scv.pp.momentsにまとめられています (内部でscv.pp.pcascv.pp.neighborsを利用)。
  • 決定論的速度推定には1次が必要ですが,確率的推定には2次モーメントが必要です.
In [14]:
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
Filtered out 20801 genes that are detected 20 counts (shared).
Normalized count data: X, spliced, unspliced.
Logarithmized X.
computing neighbors
    finished (0:00:00) --> added 
    'distances' and 'connectivities', weighted adjacency matrices (adata.obsp)
computing moments based on connectivities
    finished (0:00:00) --> added 
    'Ms' and 'Mu', moments of spliced/unspliced abundances (adata.layers)

さらなる前処理 (optional)

  • さらなる前処理(バッチ効果補正など)を使用して、不要な変動源を除去することができます。
  • 追加の前処理ステップはXにのみ影響し、スプライス/アンスプライス数には適用されないことに注意してください。

Estimate RNA velocity

  • Velocityは遺伝子発現空間におけるベクトルであり、個々の細胞の移動の「方向」と「速度」を表します。
    • Velocity = speed + direction
  • Velocitiesは、確率的(デフォルト)または決定論的 (by setting mode='deterministic')に、スプライシング動態の転写ダイナミクスをモデル化することで得られます。
  • 各遺伝子について、pre-mature (unspliced), mature (spliced)の mRNAカウントの定常状態比 (steady-state-ratio) をフィットさせ、これが一定の転写状態を表す。
  • そして、Velocityはこのsteady-state-ratioの残差として求めます。
  • 正のVelocityは、ある遺伝子が発現上昇していることを示しており、定常状態で予想されるよりも多くunspliced mRNAが存在する細胞で発生します。
  • 逆に、負の速度は遺伝子がダウンレギュレーションされていることを示している。

  • 完全な動的モデルの解は、mode='dynamical'を設定することで得られます。このモードは事前に scv.tl.recover_dynamics(adata) を実行しておく必要があります。このモデルについては次のチュートリアルで詳しく説明します。

In [15]:
scv.tl.velocity(adata)
computing velocities
    finished (0:00:01) --> added 
    'velocity', velocity vectors for each individual cell (adata.layers)

計算された速度は、カウント行列と同様に adata.layer に保存されます。

  • 遺伝子間のvelocityの組み合わせは、個々の細胞の将来の状態を推定するために使用することができます。
  • 速度を低次元の埋め込みに投影するために、細胞間遷移確率が推定されます。つまり、各速度ベクトルについて、その方向に沿ったセル遷移の可能性が高いものを見つけます。
  • 遷移確率は、潜在的な細胞間遷移と速度ベクトルとの間の余弦相関を用いて計算され、velocity graphと呼ばれる行列に格納される。
  • 結果として得られるvelocity graphは$n_{obs} \times n_{obs}$ の次元を持ち、速度ベクトルによって説明できる細胞の状態変化を要約しています(高速化のために、approx=Trueを設定することで、縮小されたPCA空間上で計算することもできます)。
In [16]:
scv.tl.velocity_graph(adata)
computing velocity graph
    finished (0:00:10) --> added 
    'velocity_graph', sparse matrix with cosine correlations (adata.uns)
  • 様々なアプリケーションのために、速度グラフを遷移行列に変換することができます。これはガウスカーネルを適用してコサイン類似度を実際の遷移確率に変換することで行います。
  • マルコフ遷移行列には scv.utils.get_transition_matrix でアクセスできます。

  • 前述したように、velocityを低次元空間に埋め込むために、scv.tl.velocity_embeddingを使用して得られた遷移確率に対する平均遷移を計算します。

  • さらに、マルコフ連鎖に沿って細胞をその起源と潜在的な運命までトレースすることができ、それにより、軌跡の根と終点を得ることができます。これは scv.tl.terminal_states で参照できます。

Project the velocities

  • 最後に、速度は任意の低次元軸に投影され、以下のいずれかの方法で可視化されます。

    • on cellular level with scv.pl.velocity_embedding,
    • as gridlines with scv.pl.velocity_embedding_grid,
    • or as streamlines with scv.pl.velocity_embedding_stream.
  • デモデータには既に計算済みの UMAP 埋め込みと注釈付きクラスタが含まれています。

  • 独自のデータに適用する場合、これらは scv.tl.umapscv.tl.louvain で取得できます。詳細については、scanpyのチュートリアルを参照してください。
  • さらに、すべてのプロット関数はデフォルトでbasis='umap'color='cluster'を使用するように設定されています。
In [19]:
scv.pl.velocity_embedding_stream(adata, basis='umap')
computing velocity embedding
    finished (0:00:00) --> added
    'velocity_umap', embedded velocity vectors (adata.obsm)
  • ストリームラインとして表示されるvelocity ベクトルフィールドは、発生過程の詳細な洞察をもたらします。
    • 管腔細胞(ductal cells) と内分泌前駆細胞 (endocrine progenitors) の周期的集団を正確に描き出します。
    • さらに、細胞の系統コミットメント、細胞周期の出口、および内分泌細胞の分化の細胞の状態を表す。

cell level velocity

  • 速度ベクトル場の最も詳細な分解能は、単一細胞レベルで得られるもので、各矢印は個々の細胞の移動方向と速度を示している。
  • これにより、例えば、Ngn3細胞の初期の内分泌関与(黄色)や、近末端のα細胞(青)と一過性のβ細胞(緑)の間の明確な違いが明らかになりました。
In [20]:
scv.pl.velocity_embedding(adata, arrow_length=3, arrow_size=2, dpi=120)

Interprete the velocities

これはおそらく最も重要な部分です。つまり、生物学的な結論を投影されたvelocityに限定するのではなく、推定された方向が特定の遺伝子にどのように支持されているかを理解するために、相図 (phase portraits)を介して個々の遺伝子のダイナミクスを調べるべきです。

spliced vs. unspliced phase portrait をどのように解釈するかについては、こちらのgif を参照してください。

  • 遺伝子の活性は転写制御によって調整されている。
  • 特定の遺伝子に対して転写が誘導されると、(新たに転写された)precursor unspliced mRNAsが増加し、逆に転写が抑制されたり、転写が行われなかったりすると、unspliced mRNAsが減少します。
  • スプライスされたmRNAは、未転写mRNAから産生され、時間差を持って同じ傾向をたどる。
  • 時間は、 hidden/latent 変数である。
  • そのため、実際に測定されているもの、すなわち、相図に表示されているspliced and unspliced mRNAsからダイナミクスを推測する必要があります。

さて、いくつかのマーカー遺伝子の相図を scv.pl.velocity(adata, gene_names) または scv.pl.scatter(adata, gene_names) で可視化してみましょう。

In [21]:
scv.pl.velocity(adata, ['Cpe',  'Gnao1', 'Ins2', 'Adk'], ncols=2)
  • 黒線は推定された ‘steady-state’ の比率、すなわち、一定の転写状態における unspliced to spliced mRNA abundanceの比率に対応しています。
  • 特定の遺伝子のRNA velocity は、残差、すなわち観測がその定常状態線 (黒線)からどれだけ逸脱しているかとして決定されます。
  • 正の速度は、ある遺伝子が発現上昇していることを示しており、定常状態で予想されるよりもunspliced mRNAが豊富な細胞で発生します。
  • 逆に、負の速度は遺伝子が発現抑制されていることを示している。

  • 例えば、Cpeは発現上昇したNgn3(黄色)からPre-endocrine(オレンジ)からβ細胞(緑)への方向性を説明し、Adkは発現抑制されたDuctal(濃い緑)からNgn3(黄色)から残りの内分泌細胞への方向性を説明している。

In [22]:
scv.pl.scatter(adata, 'Cpe',
               color=['clusters', 'velocity'],
               add_outline='Ngn3 high EP, Pre-endocrine, Beta')

Identify important genes

私たちは、結果として得られるベクターフィールドと推定される系統を説明するのに役立つ遺伝子を特定する体系的な方法を必要としています。そのために、どの遺伝子がクラスター特異的な微分速度発現を持ち、残りの集団と比較して有意に高いか低いかを検定するモジュールを提供します。

  • scv.tl.rank_velocity_genes モジュールは differential velocity t-test を実行し、各クラスタの遺伝子ランキングを出力します。
  • 遺伝子候補選択のテストを制限するための閾値を設定することができます (例: min_corr)。
In [23]:
scv.tl.rank_velocity_genes(adata, groupby='clusters', min_corr=.3)

df = scv.DataFrame(adata.uns['rank_velocity_genes']['names'])
df.head()
ranking velocity genes
    finished (0:00:02) --> added 
    'rank_velocity_genes', sorted scores by group ids (adata.uns) 
    'spearmans_score', spearmans correlation scores (adata.var)
Out[23]:
Ductal Ngn3 low EP Ngn3 high EP Pre-endocrine Beta Alpha Delta Epsilon
0 Notch2 Ptpn3 Pde1c Pam Pax6 Zcchc16 Zdbf2 Tmcc3
1 Sox5 Hacd1 Ptprs Sdk1 Unc5c Nlgn1 Spock3 Heg1
2 Krt19 Hspa8 Pclo Baiap3 Nnat Nell1 Akr1c19 Gpr179
3 Hspa8 Gm8113 Rap1gap2 Abcc8 Tmem108 Prune2 Ptprt Ica1
4 Ano6 Kcnq1 Ttyh2 Gnas Ptprt Ksr2 Snap25 Ncoa7
In [24]:
kwargs = dict(frameon=False, size=10, linewidth=1.5,
              add_outline='Ngn3 high EP, Pre-endocrine, Beta')

scv.pl.scatter(adata, df['Ngn3 high EP'][:5], ylabel='Ngn3 high EP', **kwargs)
scv.pl.scatter(adata, df['Pre-endocrine'][:5], ylabel='Pre-endocrine', **kwargs)

The genes Ptprs, Pclo, Pam, Abcc8, Gnas, for instance, support the directionality from Ngn3 high EP (yellow) to Pre-endocrine (orange) to Beta (green).

Velocities in cycling progenitors

RNA速度によって検出される細胞周期は、細胞周期スコア(位相マーカー遺伝子の平均発現量を標準化したスコア)によって生物学的に肯定されます。

In [26]:
scv.tl.score_genes_cell_cycle(adata)
scv.pl.scatter(adata, color_gradients=['S_score', 'G2M_score'], smooth=True, perc=[5, 95])
calculating cell cycle phase
-->     'S_score' and 'G2M_score', scores of cell cycle phases (adata.obs)

周期性のある管状細胞 (Ductal cells)については、SとG2Mの位相マーカーを介してスクリーニングすることができます。上のモジュールはまた、位相マーカー遺伝子のランク付けやソートに使用できるスピアマン相関スコアを計算しており、それらのphase portraitsを表示するために使用することができます。

In [27]:
s_genes, g2m_genes = scv.utils.get_phase_marker_genes(adata)
s_genes = scv.get_df(adata[:, s_genes], 'spearmans_score', sort_values=True).index
g2m_genes = scv.get_df(adata[:, g2m_genes], 'spearmans_score', sort_values=True).index

kwargs = dict(frameon=False, ylabel='cell cycle genes')
scv.pl.scatter(adata, list(s_genes[:2]) + list(g2m_genes[:3]), **kwargs)

特にHellsとTop2aは、循環する前駆体のベクトル場を説明するのに適している。Top2aは、G2M段階で実際にピークを迎える直前に高い速度が割り当てられる。そこでは、負の速度は、その後のダウンレギュレーションの直後に完全に一致する。

In [31]:
scv.pl.velocity(adata, ['Hells', 'Top2a'], ncols=2, add_outline=True)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-31-b78385210b7b> in <module>
----> 1 scv.pl.velocity(adata, ['Hells', 'Top2a'], ncols=2, add_outline=True)

/opt/conda/lib/python3.7/site-packages/scvelo/plotting/velocity.py in velocity(adata, var_names, basis, vkey, mode, fits, layers, color, color_map, colorbar, perc, alpha, size, groupby, groups, legend_loc, legend_fontsize, use_raw, fontsize, figsize, dpi, show, save, ax, ncols, **kwargs)
    122         scatter(adata, basis=var, color=color, colorbar=colorbar, color_map=cmap, perc=perc, frameon=True, title=var,
    123                 size=size, use_raw=use_raw, alpha=alpha, fontsize=fontsize, vkey=fits, legend_fontsize=legend_fontsize,
--> 124                 show=False, ax=ax, save=False, legend_loc_lines='none' if v < len(var_names)-1 else legend_loc, **kwargs)
    125 
    126         # velocity and expression plots

/opt/conda/lib/python3.7/site-packages/scvelo/plotting/scatter.py in scatter(adata, basis, x, y, vkey, color, use_raw, layer, color_map, colorbar, palette, size, alpha, linewidth, linecolor, perc, groups, sort_order, components, projection, legend_loc, legend_loc_lines, legend_fontsize, legend_fontweight, xlabel, ylabel, title, fontsize, figsize, xlim, ylim, add_density, add_assignments, add_linfit, add_polyfit, add_rug, add_text, add_text_pos, add_outline, outline_width, outline_color, n_convolve, smooth, rescale_color, color_gradients, dpi, frameon, zorder, ncols, nrows, wspace, hspace, show, save, ax, **kwargs)
    375                     if add_outline in adata.var.keys() and basis in adata.var_names:
    376                         add_outline = str(adata[:, basis].var[add_outline][0])
--> 377                 idx = groups_to_bool(adata, add_outline, color)
    378                 if idx is not None and np.sum(idx) > 0:  # if anything to be outlined
    379                     zorder = 2 if zorder is None else zorder + 2

/opt/conda/lib/python3.7/site-packages/scvelo/plotting/utils.py in groups_to_bool(adata, groups, groupby)
    168 
    169 def groups_to_bool(adata, groups, groupby=None):
--> 170     groups, groupby = get_groups(adata, groups, groupby)
    171     if isinstance(groups, (list, tuple, np.ndarray, np.record)):
    172         if groupby is not None and isinstance(groups[0], str):

TypeError: cannot unpack non-iterable NoneType object

Speed and coherence

さらに2つの有用な統計量があります。

  • 微分の速度や速度は速度ベクトルの長さで与えられます。
  • ベクトル場のコヒーレンス(すなわち、速度ベクトルが隣接する速度とどのように相関するか)は、信頼性の尺度となります。
In [32]:
scv.tl.velocity_confidence(adata)
keys = 'velocity_length', 'velocity_confidence'
scv.pl.scatter(adata, c=keys, cmap='coolwarm', perc=[5, 95])
--> added 'velocity_length' (adata.obs)
--> added 'velocity_confidence' (adata.obs)
--> added 'velocity_confidence_transition' (adata.obs)

これらの結果から、細胞の分化のペースが遅い場合と早い場合、またその方向が不明な場合があることが明らかになった。

クラスターレベルでは、細胞周期の終了後(Ngn3低EP)には分化が実質的に加速され、ベータ細胞の産生中にはペースを維持し、アルファ細胞の産生中にはペースが低下することがわかった。

In [33]:
df = adata.obs.groupby('clusters')[keys].mean().T
df.style.background_gradient(cmap='coolwarm', axis=1)
Out[33]:
clusters Ductal Ngn3 low EP Ngn3 high EP Pre-endocrine Beta Alpha Delta Epsilon
velocity_length 5.707904 6.227023 13.257741 13.429206 14.219882 10.491913 7.623143 10.477606
velocity_confidence 0.775042 0.756474 0.896944 0.830438 0.725973 0.752715 0.691188 0.784467

Velocity graph and pseudotime

  • Velocityグラフを可視化して、推定された全ての細胞間接続/遷移を描写することができます。
  • 閾値を設定することで、高確率の遷移に限定することができる。
  • 例えば、このグラフは、初期と後期の前内分泌細胞からのイプシロン細胞生産の2つのフェーズを示しています。
In [34]:
scv.pl.velocity_graph(adata, threshold=.1)

さらに、このグラフは、指定された細胞の子孫/先祖を描くために使用することができます。ここでは、前内分泌細胞は、その潜在的な運命までトレースされています。

In [35]:
x, y = scv.utils.get_cell_transitions(adata, basis='umap', starting_cell=70)
ax = scv.pl.velocity_graph(adata, c='lightgrey', edge_width=.05, show=False)
ax = scv.pl.scatter(adata, x=x, y=y, s=120, c='ascending', cmap='gnuplot', ax=ax)
  • 最後に、速度グラフに基づいて、 velocity pseudotimeを計算することができる。
  • グラフからroot細胞からの分布を推定した後、root細胞からグラフに沿って歩いた後、その細胞に到達するまでにかかる平均ステップ数を測定することで計算されます。

  • diffusion pseudotimeとは逆に、root細胞を暗黙的に推定し、similarity-based diffusion kernelの代わりに有向速度グラフに基づいています。

In [36]:
scv.tl.velocity_pseudotime(adata)
scv.pl.scatter(adata, color='velocity_pseudotime', cmap='gnuplot')
computing terminal states
    identified 2 regions of root cells and 1 region of end points 
    finished (0:00:00) --> added
    'root_cells', root cells of Markov diffusion process (adata.obs)
    'end_points', end points of Markov diffusion process (adata.obs)

PAGA velocity graph

  • PAGAグラフ抽象化は、軌跡推論のための最高性能の手法としてベンチマーキングされています。
  • PAGAは、2つのクラスタ間の接続性に対応した重み付きエッジを用いて、データのトポロジーをグラフ化したようなマップを提供する。
  • ここでは、PAGAを速度推定の方向性によって拡張している。
In [39]:
scv.tl.paga(adata, groups='clusters')
df = scv.get_df(adata, 'paga/transitions_confidence', precision=2).T
df.style.background_gradient(cmap='Blues').format('{:.2g}')
running PAGA
    finished (0:00:02) --> added
    'paga/transitions_confidence', connectivities adjacency (adata.uns)
    'paga/connectivities', connectivities adjacency (adata.uns)
    'paga/connectivities_tree', connectivities subtree (adata.uns)
Out[39]:
Ductal Ngn3 low EP Ngn3 high EP Pre-endocrine Beta Alpha Delta Epsilon
Ductal 0 0.15 0 0 0 0 0 0
Ngn3 low EP 0 0 0.24 0 0 0 0 0
Ngn3 high EP 0 0 0 0.23 0 0 0 0
Pre-endocrine 0 0 0 0 0.5 0.12 0.21 0.12
Beta 0 0 0 0 0 0 0 0
Alpha 0 0 0 0 0 0 0 0
Delta 0 0 0 0 0 0 0 0
Epsilon 0 0 0 0 0 0 0 0

これは、左/行から右/列へと読み取ることができ、例えば、DuctalからNgn3 low EPへの信頼性ある遷移を割り当てることができます。

この表は、UMAP上に埋め込まれた有向グラフでまとめることができます。

In [40]:
scv.pl.paga(adata, basis='umap',
            size=50,
            alpha=.1,
            min_edge_width=2,
            node_size_scale=1.5)
In [ ]: